load required packages

library(qiime2R)
library(tidyverse)
library(ggplot2)
library(lubridate)
library(xts)
library(ggpubr)
library(dplyr)
library(data.table)
library(vegan)
#library(reshape2)
#install.packages("janitor")
#library(janitor)
library(viridis)

show versions:

sessionInfo()
## R version 4.3.2 (2023-10-31)
## Platform: x86_64-apple-darwin20 (64-bit)
## Running under: macOS Monterey 12.2.1
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRblas.0.dylib 
## LAPACK: /Library/Frameworks/R.framework/Versions/4.3-x86_64/Resources/lib/libRlapack.dylib;  LAPACK version 3.11.0
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## time zone: Europe/London
## tzcode source: internal
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] viridis_0.6.4     viridisLite_0.4.2 vegan_2.6-4       lattice_0.21-9   
##  [5] permute_0.9-7     data.table_1.15.0 ggpubr_0.6.0      xts_0.13.2       
##  [9] zoo_1.8-12        lubridate_1.9.3   forcats_1.0.0     stringr_1.5.1    
## [13] dplyr_1.1.4       purrr_1.0.2       readr_2.1.4       tidyr_1.3.1      
## [17] tibble_3.2.1      ggplot2_3.4.4     tidyverse_2.0.0   qiime2R_0.99.6   
## 
## loaded via a namespace (and not attached):
##  [1] bitops_1.0-7            gridExtra_2.3           rlang_1.1.3            
##  [4] magrittr_2.0.3          ade4_1.7-22             compiler_4.3.2         
##  [7] mgcv_1.9-0              vctrs_0.6.5             reshape2_1.4.4         
## [10] pkgconfig_2.0.3         crayon_1.5.2            fastmap_1.1.1          
## [13] backports_1.4.1         XVector_0.42.0          utf8_1.2.4             
## [16] rmarkdown_2.25          tzdb_0.4.0              xfun_0.42              
## [19] zlibbioc_1.48.0         cachem_1.0.8            GenomeInfoDb_1.38.4    
## [22] jsonlite_1.8.8          biomformat_1.30.0       rhdf5filters_1.14.1    
## [25] Rhdf5lib_1.24.1         broom_1.0.5             parallel_4.3.2         
## [28] cluster_2.1.4           R6_2.5.1                bslib_0.6.1            
## [31] stringi_1.8.3           car_3.1-2               zCompositions_1.5      
## [34] rpart_4.1.21            jquerylib_0.1.4         Rcpp_1.0.12            
## [37] iterators_1.0.14        knitr_1.45              base64enc_0.1-3        
## [40] IRanges_2.36.0          Matrix_1.6-1.1          splines_4.3.2          
## [43] nnet_7.3-19             igraph_1.6.0            timechange_0.2.0       
## [46] tidyselect_1.2.0        abind_1.4-5             rstudioapi_0.15.0      
## [49] yaml_2.3.8              codetools_0.2-19        plyr_1.8.9             
## [52] Biobase_2.62.0          withr_3.0.0             evaluate_0.23          
## [55] foreign_0.8-85          survival_3.5-7          Biostrings_2.70.1      
## [58] pillar_1.9.0            phyloseq_1.46.0         carData_3.0-5          
## [61] checkmate_2.3.1         DT_0.31                 foreach_1.5.2          
## [64] stats4_4.3.2            NADA_1.6-1.1            generics_0.1.3         
## [67] RCurl_1.98-1.13         truncnorm_1.0-9         S4Vectors_0.40.2       
## [70] hms_1.1.3               munsell_0.5.0           scales_1.3.0           
## [73] glue_1.7.0              Hmisc_5.1-1             tools_4.3.2            
## [76] ggsignif_0.6.4          rhdf5_2.46.1            grid_4.3.2             
## [79] ape_5.7-1               colorspace_2.1-0        nlme_3.1-163           
## [82] GenomeInfoDbData_1.2.11 htmlTable_2.4.2         Formula_1.2-5          
## [85] cli_3.6.2               fansi_1.0.6             gtable_0.3.4           
## [88] rstatix_0.7.2           sass_0.4.8              digest_0.6.34          
## [91] BiocGenerics_0.48.1     htmlwidgets_1.6.4       htmltools_0.5.7        
## [94] multtest_2.58.0         lifecycle_1.0.4         MASS_7.3-60

Eukaryotes

18S metabarcoding data analysis ## Beta Diversity

setwd("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/")

SVs <- read_qza("rarefied_table.qza")$data
metadata <- read_q2metadata("metadata-rarefied.txt")
taxonomy <- read_qza("taxonomy.qza")$data %>% parse_taxonomy()
taxonomy$Kingdom <- gsub("d__", "",taxonomy$Kingdom)

metadata <- data.frame(metadata[,-1], row.names=metadata[,1])

taxasums<-summarize_taxa(SVs, taxonomy)$Species

# now I'm going to create betadiversity plots and look at the what environmental
# variables may be influencing beta diversity
mtaxasums <- t(as.matrix(taxasums))
set.seed(123)
nmds <- metaMDS(mtaxasums, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1898539 
## Run 1 stress 0.1947476 
## Run 2 stress 0.1918644 
## Run 3 stress 0.1998882 
## Run 4 stress 0.194419 
## Run 5 stress 0.1952102 
## Run 6 stress 0.1916873 
## Run 7 stress 0.200839 
## Run 8 stress 0.2079919 
## Run 9 stress 0.1994643 
## Run 10 stress 0.3982293 
## Run 11 stress 0.1940355 
## Run 12 stress 0.1943068 
## Run 13 stress 0.2128468 
## Run 14 stress 0.18729 
## ... New best solution
## ... Procrustes: rmse 0.05572826  max resid 0.1936431 
## Run 15 stress 0.2062908 
## Run 16 stress 0.1991352 
## Run 17 stress 0.2159338 
## Run 18 stress 0.189858 
## Run 19 stress 0.1876837 
## ... Procrustes: rmse 0.02183372  max resid 0.09050004 
## Run 20 stress 0.1940356 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
en <- envfit(nmds, metadata, permutations = 999, na.rm = TRUE)

#extract NMDS scores (x and y coordinates) for sites from newer versions of vegan package
data.scores <- as.data.frame(scores(nmds)$sites)
data.score <- merge(data.scores, metadata, by=0) # adding metadata
#extracting vectors
en_coord_cont = as.data.frame(scores(en, "vectors")) * ordiArrowMul(en)
en_coord_cat = as.data.frame(scores(en, "factors")) * ordiArrowMul(en)

Mantel Tests

# To work out which environmental vectors to include I'm going to do some mantel tests
# first create the distance matrices then do the mantel test
dist.abundance <- vegdist(t(taxasums), method = "bray")
dist.time <- dist(metadata$Time.Since.Start, method = "euclidean")
abund_time = mantel(dist.abundance, dist.time, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_time
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.time, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.1511 
##       Significance: 0.005 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0700 0.0931 0.1112 0.1332 
## Permutation: free
## Number of permutations: 9999
dist.day.night <- dist(metadata$nDay_Night, method = "euclidean")
abund_day.night = mantel(dist.abundance, dist.day.night, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_day.night
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.day.night, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: -0.03779 
##       Significance: 0.7517 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0688 0.0903 0.1081 0.1353 
## Permutation: free
## Number of permutations: 9999
dist.moon <- dist(metadata$nMoon_Moonset, method = "euclidean")
abund_moon = mantel(dist.abundance, dist.moon, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_moon
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.moon, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: -0.09146 
##       Significance: 0.9681 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0727 0.0956 0.1166 0.1413 
## Permutation: free
## Number of permutations: 9999
dist.tide <- dist(metadata$UKHO_Tide_Height, method = "euclidean")
abund_tide = mantel(dist.abundance, dist.tide, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_tide
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.tide, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: -0.01736 
##       Significance: 0.6145 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0638 0.0853 0.1053 0.1244 
## Permutation: free
## Number of permutations: 9999
dist.caltemp <- dist(metadata$Calista_Temp_oC, method = "euclidean")
abund_caltemp = mantel(dist.abundance, dist.caltemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_caltemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.caltemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.2309 
##       Significance: 0.0026 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0976 0.1295 0.1567 0.1846 
## Permutation: free
## Number of permutations: 9999
dist.airtemp <- dist(metadata$SotonMet_Air_Temp_oC, method = "euclidean")
abund_airtemp = mantel(dist.abundance, dist.airtemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_airtemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.airtemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.02534 
##       Significance: 0.3255 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0858 0.1127 0.1373 0.1621 
## Permutation: free
## Number of permutations: 9999
dist.wind <- dist(metadata$SotonMet_Wind_Speed_knots, method = "euclidean")
abund_wind = mantel(dist.abundance, dist.wind, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_wind
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.wind, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.132 
##       Significance: 0.0801 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.118 0.156 0.187 0.224 
## Permutation: free
## Number of permutations: 9999
dist.winddirection <- dist(metadata$SotonMet_Wind_Direction_degrees, method = "euclidean")
abund_winddirection = mantel(dist.abundance, dist.winddirection, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_winddirection
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.winddirection, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.101 
##       Significance: 0.1405 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.121 0.158 0.193 0.235 
## Permutation: free
## Number of permutations: 9999
dist.gust <- dist(metadata$SotonMet_Max_Gust_knots, method = "euclidean")
abund_gust = mantel(dist.abundance, dist.gust, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_gust
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.gust, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.07845 
##       Significance: 0.1727 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.109 0.146 0.179 0.214 
## Permutation: free
## Number of permutations: 9999
dist.pressure <- dist(metadata$SotonMet_Barometric_Pressure_millibars, method = "euclidean")
abund_pressure = mantel(dist.abundance, dist.pressure, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_pressure
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.pressure, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.1275 
##       Significance: 0.0933 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.123 0.163 0.199 0.239 
## Permutation: free
## Number of permutations: 9999
dist.UoStemp <- dist(metadata$UoS_Temp_oC, method = "euclidean")
abund_UoStemp = mantel(dist.abundance, dist.UoStemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_UoStemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.UoStemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.1409 
##       Significance: 0.0455 
## 
## Upper quantiles of permutations (null model):
##   90%   95% 97.5%   99% 
## 0.104 0.136 0.167 0.195 
## Permutation: free
## Number of permutations: 9999

nMDS

#removing non significant variables and time since start
drops <- c("Conc_NOC_ngµl",
           "Conc_ESS_ngµl",
           "Conc_PCR2",
           "Conc_diff_ngµl",
           "nDay_Night",
           "nMoon_Moonset",
           "SampleNo",
           "UKHO_Tide_Height",
           "SotonMet_Air_Temp_oC",
           "SotonMet_Wind_Speed_knots",
           "SotonMet_Wind_Direction_degrees",
           "SotonMet_Max_Gust_knots",
           "SotonMet_Barometric_Pressure_millibars",
           "UoS_Temp_oC", 
           "Time.Since.Start")
en_coord_cont <- en_coord_cont[!(row.names(en_coord_cont) %in% drops),]

nMDS18S = ggplot(data = data.score, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(data = data.score, aes(color=`Time.Since.Start`/60), size = 3, alpha = 0.8) + 
  scale_fill_viridis(limits = c(0, 94))  +
  scale_colour_viridis(limits = c(0, 94), 
                        breaks = c(0, 24, 48, 72),
                        labels = c(0, 24, 48, 72)) +
  geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
               data = en_coord_cont, linewidth =1, alpha = 0.5, colour = "grey30") +
  geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), colour = "grey30", 
            fontface = "bold", label = "Temp", parse = TRUE, position = position_nudge(x = -0.05)) + 
  theme(axis.title = element_text(size = 10, face = "bold", colour = "grey30"), 
        panel.background = element_blank(), panel.border = element_rect(fill = NA, colour = "grey30"), 
        axis.ticks = element_blank(), axis.text = element_blank(), legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "grey30"), 
        legend.text = element_text(size = 9, colour = "grey30")) +
  theme_q2r() +
  labs(title="Eukaryotes", colour = "Hours Since Start")
nMDS18S

## Alpha Diversity Time Series

shannon<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/core-metrics-results/shannon_vector.qza")
shannon<-shannon$data %>% rownames_to_column("SampleID") # this moves the sample names to a new column that matches the metadata and allows them to be merged

observedASVs<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/core-metrics-results/observed_features_vector.qza")
observedASVs<-observedASVs$data %>% rownames_to_column("SampleID")

evenness<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/core-metrics-results/evenness_vector.qza")
evenness<-evenness$data %>% rownames_to_column("SampleID")

faithPD<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/core-metrics-results/faith_pd_vector.qza")
faithPD<-faithPD$data %>% rownames_to_column("SampleID")

metadata<- metadata %>% rownames_to_column("SampleID") # I think this just assigned 1,2,3,4
# instead I need to change the Sample column to SampleID
## first remove the new sample ID column
# metadata = subset(metadata, select = -SampleID)
# then rename Sample to SampleID
# names(metadata)[names(metadata) == 'sample'] <- 'SampleID'

#assign alpha diversity metrics to the metadata
metadata<-
  metadata %>% 
  left_join(shannon) %>%
  left_join(observedASVs) %>%
  left_join(evenness) %>%
  left_join(faithPD)
## Joining with `by = join_by(SampleID)`
## Joining with `by = join_by(SampleID)`
## Joining with `by = join_by(SampleID)`
## Joining with `by = join_by(SampleID)`
metadata$`Start_Date_Time`<-dmy_hm(metadata$`Start_Date_Time`) # getting the dates into ISO format

ShannonPlot <- ggplot(metadata, aes(x = `Start_Date_Time`, y=`shannon_entropy`)) +
  geom_line() +
  geom_point() +
  ylim(2,7.5) +
  theme_q2r() +
  xlab("Time") +
  ylab("Diversity")
ShannonPlot 

EvennessPlot <- ggplot(metadata, aes(x = `Start_Date_Time`, y=`pielou_evenness`)) +
  geom_line() +
  geom_point() +
  ylim(0.25,0.85) +
  theme_q2r() +
  xlab("Time") +
  ylab("Evenness")
EvennessPlot 

ObservedPlot <- ggplot(metadata, aes(x = `Start_Date_Time`, y=`observed_features`)) +
  geom_line() +
  geom_point() +
  theme_q2r() +
  xlab("Time") +
  ylab("ASVs") +
  labs(title = "Eukaryotes")
ObservedPlot

FaithPlot <- ggplot(metadata, aes(x = `Start_Date_Time`, y=`faith_pd`)) +
  geom_line() +
  geom_point() +
  theme_q2r() +
  xlab("Time") +
  ylab("Faith PD")
FaithPlot 

Alpha18STime<- ggarrange(ObservedPlot, ShannonPlot, EvennessPlot, FaithPlot,
          ncol = 1, 
          align = "hv")
Alpha18STime

## Stacked Taxa Time Series ### Species

setwd("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/18S_dada2redo/")
metadata <- read_q2metadata("metadata-rarefied.txt")
metadata <- data.frame(metadata[,-1], row.names=metadata[,1])
taxasums <- cbind(taxasums, total = rowSums(taxasums)) #add total col for each species
taxasums <- taxasums[order(-taxasums$total),] #order by total column
write.csv(taxasums, "SpeciesList.csv")

# seperate the top 10 taxa
top10species <- taxasums %>%
  arrange(desc(total)) %>%
  slice(1:10)
#remove the top 10 species 
remainderSpecies <- taxasums[-(1:10),]
remainderSpecies <- rbind(remainderSpecies, Remainder = colSums(remainderSpecies)) #add total row for each sample
# rejoin the tables
taxasums10 <- rbind(top10species, remainderSpecies)
# take the top 11 which will now include the remainder total
taxasums10 <- taxasums10 %>%
  arrange(desc(total)) %>%
  slice(1:11)
#remove the total column
taxasums10$total <- NULL
#rownames to first column
taxasums10<- tibble::rownames_to_column(taxasums10)
taxasums10<- data.frame(t(taxasums10)) #transpose
names(taxasums10)<- taxasums10[1,] # make the taxa the rownames
taxasums10<- taxasums10[-1,] # remove the first duplicate first row

# add Date_Time to taxasums10
t.metadata <- metadata[,c("sample","Start_Date_Time")]
t.metadata$Start_Date_Time <- as.character(t.metadata$Start_Date_Time) #convert posixct back to character



#t.metadata <- row_to_names(t.metadata, row_number = 1, remove_row = TRUE) #make the sample names col names
taxasums101 <- cbind(t.metadata, taxasums10) #combine taxa and times
taxasums101$Start_Date_Time<- dmy_hm(taxasums101$Start_Date_Time)
write.csv(taxasums101, "species10.csv")

taxasums101$sample<- NULL # remove sample column

species<- taxasums101 %>%
  pivot_longer(!Start_Date_Time,, names_to = "Taxa", values_to = "Sequence count")

species$`Sequence count`<- as.numeric(species$`Sequence count`)

#plot stack
species18SStack <- species %>%
  ggplot(aes(fill= Taxa,
             x=Start_Date_Time,
             y=`Sequence count`,
             colour = Taxa)) +
  geom_area(position="stack", 
            stat="identity") +
  geom_line(position = 'stack', linewidth = 1, colour = "black" ) +
  scale_colour_viridis(discrete = TRUE, direction = -1) +
  scale_fill_viridis(discrete = TRUE, direction = -1) +
  theme_q2r() +
  theme(legend.position = "bottom", legend.box = "vertical", legend.justification = "left", legend.box.just = "left", legend.direction = "vertical") +
  labs(title = "Eukaryotes", x="Time")
  

species18SStack

### Genus

Prokaryotes

Beta Diversity

setwd("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/")

SVs <- read_qza("rarefied_table.qza")$data
metadata <- read_q2metadata("metadata-no-blanks.txt")
taxonomy <- read_qza("taxonomy.qza")$data %>% parse_taxonomy()
taxonomy$Kingdom <- gsub("d__", "",taxonomy$Kingdom)

metadata <- data.frame(metadata[,-1], row.names=metadata[,1])

taxasums<-summarize_taxa(SVs, taxonomy)$Species

# now I'm going to create betadiversity plots and look at the what environmental
# variables may be influencing beta diversity
mtaxasums <- t(as.matrix(taxasums))
set.seed(123)
nmds <- metaMDS(mtaxasums, distance = "bray")
## Square root transformation
## Wisconsin double standardization
## Run 0 stress 0.1929938 
## Run 1 stress 0.2104758 
## Run 2 stress 0.2122065 
## Run 3 stress 0.2017629 
## Run 4 stress 0.203068 
## Run 5 stress 0.2175578 
## Run 6 stress 0.207389 
## Run 7 stress 0.398229 
## Run 8 stress 0.1902876 
## ... New best solution
## ... Procrustes: rmse 0.07528947  max resid 0.3149976 
## Run 9 stress 0.1891515 
## ... New best solution
## ... Procrustes: rmse 0.03557354  max resid 0.1759746 
## Run 10 stress 0.1891534 
## ... Procrustes: rmse 0.003089313  max resid 0.01691136 
## Run 11 stress 0.1952534 
## Run 12 stress 0.1907372 
## Run 13 stress 0.1952534 
## Run 14 stress 0.1952537 
## Run 15 stress 0.2033445 
## Run 16 stress 0.1888986 
## ... New best solution
## ... Procrustes: rmse 0.01907387  max resid 0.08086317 
## Run 17 stress 0.1909827 
## Run 18 stress 0.1918261 
## Run 19 stress 0.1937858 
## Run 20 stress 0.2114092 
## *** Best solution was not repeated -- monoMDS stopping criteria:
##     20: stress ratio > sratmax
en <- envfit(nmds, metadata, permutations = 999, na.rm = TRUE)

# base r plots are a mess so going to use ggplot2
#extract NMDS scores (x and y coordinates) for sites from newer versions of vegan package
data.scores <- as.data.frame(scores(nmds)$sites)
data.score <- merge(data.scores, metadata, by=0) # adding metadata
#extracting vectors
en_coord_cont = as.data.frame(scores(en, "vectors")) * ordiArrowMul(en)
en_coord_cat = as.data.frame(scores(en, "factors")) * ordiArrowMul(en)

Mantel Tests

# To work out which environmental vectors to include I'm going to do some mantel tests
# first create the distance matrices then do the mantel test
dist.abundance <- vegdist(t(taxasums), method = "bray")

dist.time <- dist(metadata$Time.Since.Start, method = "euclidean")
abund_time = mantel(dist.abundance, dist.time, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_time
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.time, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.5896 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0600 0.0809 0.1001 0.1265 
## Permutation: free
## Number of permutations: 9999
dist.day.night <- dist(metadata$nDay_Night, method = "euclidean")
abund_day.night = mantel(dist.abundance, dist.day.night, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_day.night
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.day.night, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: -0.06712 
##       Significance: 0.9412 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0605 0.0815 0.1013 0.1267 
## Permutation: free
## Number of permutations: 9999
dist.moon <- dist(metadata$nMoon_Moonset, method = "euclidean")
abund_moon = mantel(dist.abundance, dist.moon, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_moon
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.moon, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: -0.02434 
##       Significance: 0.6761 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0622 0.0850 0.1042 0.1305 
## Permutation: free
## Number of permutations: 9999
dist.tide <- dist(metadata$UKHO_Tide_Height, method = "euclidean")
abund_tide = mantel(dist.abundance, dist.tide, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_tide
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.tide, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.01872 
##       Significance: 0.2957 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0553 0.0755 0.0948 0.1187 
## Permutation: free
## Number of permutations: 9999
dist.caltemp <- dist(metadata$Calista_Temp_oC, method = "euclidean")
abund_caltemp = mantel(dist.abundance, dist.caltemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_caltemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.caltemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.3094 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0783 0.1022 0.1257 0.1508 
## Permutation: free
## Number of permutations: 9999
dist.airtemp <- dist(metadata$SotonMet_Air_Temp_oC, method = "euclidean")
abund_airtemp = mantel(dist.abundance, dist.airtemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_airtemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.airtemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.02901 
##       Significance: 0.2779 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0694 0.0921 0.1152 0.1457 
## Permutation: free
## Number of permutations: 9999
dist.wind <- dist(metadata$SotonMet_Wind_Speed_knots, method = "euclidean")
abund_wind = mantel(dist.abundance, dist.wind, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_wind
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.wind, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.06955 
##       Significance: 0.1595 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0943 0.1225 0.1476 0.1779 
## Permutation: free
## Number of permutations: 9999
dist.winddirection <- dist(metadata$SotonMet_Wind_Direction_degrees, method = "euclidean")
abund_winddirection = mantel(dist.abundance, dist.winddirection, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_winddirection
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.winddirection, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.4445 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0944 0.1240 0.1488 0.1825 
## Permutation: free
## Number of permutations: 9999
dist.gust <- dist(metadata$SotonMet_Max_Gust_knots, method = "euclidean")
abund_gust = mantel(dist.abundance, dist.gust, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_gust
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.gust, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.05695 
##       Significance: 0.1869 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0884 0.1174 0.1442 0.1723 
## Permutation: free
## Number of permutations: 9999
dist.pressure <- dist(metadata$SotonMet_Barometric_Pressure_millibars, method = "euclidean")
abund_pressure = mantel(dist.abundance, dist.pressure, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_pressure
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.pressure, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.2172 
##       Significance: 0.0049 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0983 0.1303 0.1606 0.1946 
## Permutation: free
## Number of permutations: 9999
dist.UoStemp <- dist(metadata$UoS_Temp_oC, method = "euclidean")
abund_UoStemp = mantel(dist.abundance, dist.UoStemp, method = "spearman", permutations = 9999, na.rm = TRUE)
abund_UoStemp
## 
## Mantel statistic based on Spearman's rank correlation rho 
## 
## Call:
## mantel(xdis = dist.abundance, ydis = dist.UoStemp, method = "spearman",      permutations = 9999, na.rm = TRUE) 
## 
## Mantel statistic r: 0.3463 
##       Significance: 1e-04 
## 
## Upper quantiles of permutations (null model):
##    90%    95%  97.5%    99% 
## 0.0833 0.1088 0.1306 0.1561 
## Permutation: free
## Number of permutations: 9999

nMDS

#removing non significant variables and time since start
drops <- c("Conc.NOC",
           "Conc.ESS",
           "Conc.PCR2",
           "Conc.diff",
           "nDay_Night",
           "nMoon_Moonset",
           "SampleNo",
           "UKHO_Tide_Height",
           "Calista_Temp_oC",
           "SotonMet_Air_Temp_oC",
           "SotonMet_Wind_Speed_knots",
           "SotonMet_Max_Gust_knots",
           "SotonMet_Wind_Direction_degrees",
           "Time.Since.Start")
en_coord_cont <- en_coord_cont[!(row.names(en_coord_cont) %in% drops),]
en_label=c("P[b]","Temp")

# now to plot and fit the environmental data
nMDS16S = ggplot(data = data.score, aes(x = NMDS1, y = NMDS2)) + 
  geom_point(data = data.score, aes(colour = `Time.Since.Start`/60), size = 3, alpha = 0.8) + 
  scale_fill_viridis(limits = c(0, 94))  +
  scale_colour_viridis(limits = c(0, 94), 
                        breaks = c(0, 24, 48, 72),
                        labels = c(0, 24, 48, 72)) +
  geom_segment(aes(x = 0, y = 0, xend = NMDS1, yend = NMDS2), 
               data = en_coord_cont, linewidth =1, alpha = 0.5, colour = "grey30") +
  geom_text(data = en_coord_cont, aes(x = NMDS1, y = NMDS2), colour = "grey30", 
            fontface = "bold", label = en_label , parse = TRUE, position = position_nudge(x = -0.05)) + 
  theme(axis.title = element_text(size = 10, face = "bold", colour = "grey30"), 
        panel.background = element_blank(), panel.border = element_rect(fill = NA, colour = "grey30"), 
        axis.ticks = element_blank(), axis.text = element_blank(), legend.key = element_blank(), 
        legend.title = element_text(size = 10, face = "bold", colour = "grey30"), 
        legend.text = element_text(size = 9, colour = "grey30")) +
  theme_q2r() +
  labs(title="Prokaryotes", colour = "Hours Since Start")
nMDS16S

## Alpha Diversity Time Series

shannon<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/core-metrics-results/shannon_vector.qza")
shannon<-shannon$data %>% rownames_to_column("Sample") # this moves the sample names to a new column that matches the metadata and allows them to be merged

observedASVs<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/core-metrics-results/observed_features_vector.qza")
observedASVs<-observedASVs$data %>% rownames_to_column("Sample")

evenness<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/core-metrics-results/evenness_vector.qza")
evenness<-evenness$data %>% rownames_to_column("Sample")

faithPD<-read_qza("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/core-metrics-results/faith_pd_vector.qza")
faithPD<-faithPD$data %>% rownames_to_column("Sample")

#assign alpha diversity metrics to the metadata
metadata<-
  metadata %>% 
  left_join(shannon) %>%
  left_join(observedASVs) %>%
  left_join(evenness) %>%
  left_join(faithPD)
## Joining with `by = join_by(Sample)`
## Joining with `by = join_by(Sample)`
## Joining with `by = join_by(Sample)`
## Joining with `by = join_by(Sample)`
metadata$`Start.Date.Time`<-dmy_hm(metadata$`Start.Date.Time`) # getting the dates into ISO format

ShannonPlot <- ggplot(metadata, aes(x = `Start.Date.Time`, y=`shannon_entropy`)) +
  geom_line() +
  geom_point() +
  ylim(2,7.5) +
  theme_q2r() +
  xlab("Time") +
  ylab("Diversity")
ShannonPlot 

EvennessPlot <- ggplot(metadata, aes(x = `Start.Date.Time`, y=`pielou_evenness`)) +
  geom_line() +
  geom_point() +
  ylim(0.25,0.85) +
  theme_q2r() +
  xlab("Time") +
  ylab("Evenness")
EvennessPlot 

ObservedPlot <- ggplot(metadata, aes(x = `Start.Date.Time`, y=`observed_features`)) +
  geom_line() +
  geom_point() +
  theme_q2r() +
  xlab("Time") +
  ylab("ASVs") +
  labs(title = "Prokaryotes")
ObservedPlot

FaithPlot <- ggplot(metadata, aes(x = `Start.Date.Time`, y=`faith_pd`)) +
  geom_line() +
  geom_point() +
  theme_q2r() +
  xlab("Time") +
  ylab("Faith PD")
FaithPlot 

Alpha16STime<- ggarrange(ObservedPlot, ShannonPlot, EvennessPlot, FaithPlot,
          ncol = 1,
          align = "hv")

Stacked Taxa Time Series

Species

# now I'm going to create a summary of the top 10 taxa and the remainder of taxa
setwd("/Users/RMS1U18/Library/CloudStorage/OneDrive-UniversityofSouthampton/PhD/Soton/Bioinformatics/16S/")
metadata <- read_q2metadata("metadata-no-blanks.txt")
metadata <- data.frame(metadata[,-1], row.names=metadata[,1])
taxasums <- cbind(taxasums, total = rowSums(taxasums)) #add total col for each species
taxasums <- taxasums[order(-taxasums$total),]
write.csv(taxasums, "SpeciesList.csv")
# seperate the top 10 taxa
top10species <- taxasums %>%
  arrange(desc(total)) %>%
  slice(1:10)
#remove the top 10 species 
remainderSpecies <- taxasums[-(1:10),]
remainderSpecies <- rbind(remainderSpecies, Remainder = colSums(remainderSpecies)) #add total row for each sample
# rejoin the tables
taxasums10 <- rbind(top10species, remainderSpecies)
# take the top 11 which will now include the remainder total
taxasums10 <- taxasums10 %>%
  arrange(desc(total)) %>%
  slice(1:11)
#remove the total column
taxasums10$total <- NULL
#rownames to first column
taxasums10<- tibble::rownames_to_column(taxasums10)
taxasums10<- data.frame(t(taxasums10)) #transpose
names(taxasums10)<- taxasums10[1,] # make the taxa the rownames
taxasums10<- taxasums10[-1,] # remove the first duplicate first row

# add Date_Time to taxasums10
t.metadata <- metadata[,c("Sample","Start.Date.Time")]
t.metadata$Start.Date.Time <- as.character(t.metadata$Start.Date.Time) #convert posixct back to character

#t.metadata <- row_to_names(t.metadata, row_number = 1, remove_row = TRUE) #make the sample names col names
taxasums101 <- cbind(t.metadata, taxasums10) #combine taxa and times
taxasums101$Start.Date.Time<- dmy_hm(taxasums101$Start.Date.Time)
write.csv(taxasums101, "species10.csv")

taxasums101$Sample<- NULL # remove sample column

species<- taxasums101 %>%
  pivot_longer(!Start.Date.Time,, names_to = "Taxa", values_to = "Sequence count")


species$`Sequence count`<- as.numeric(species$`Sequence count`)

#plot stack
species16Sstack <- species %>%
  ggplot(aes(fill= Taxa,
             x=Start.Date.Time,
             y=`Sequence count`,
             colour = Taxa)) +
  geom_area(position="stack", 
            stat="identity") +
  geom_line(position = 'stack', linewidth = 1, colour = "black" ) +
  scale_colour_viridis(discrete = TRUE, direction = -1) +
  scale_fill_viridis(discrete = TRUE, direction = -1) +
  theme_q2r() +
  theme(legend.position = "bottom", legend.box = "vertical", legend.justification = "left", legend.box.just = "left", legend.direction = "vertical") +
  labs(title = "Prokaryotes", x="Time")
  

species16Sstack

### Genus

Environmental Time Series

# set directory
setwd("~/OneDrive - University of Southampton/PhD/Soton/Soton_metadata/")
#import UoS temp data
UoStemp2019<-read.csv("UoSbuoy2019Tempdata_Cut.csv",
           stringsAsFactors = FALSE)

#import calista temp data
caltemp<-read.csv("calista_temp_dock_only.csv",
           stringsAsFactors = FALSE)

#import UKHO time an moon data
UKHO<-read.csv("UKHO_tidesunmoondata.csv",
                  stringsAsFactors = FALSE)

#import sotonmet dat inc. tide height (use stringAsFactors= FALSE to stop the data_time being read as a factor)
sotmet11<-read.csv("Sotonmet11Nov2019.csv",
                   stringsAsFactors = FALSE)
sotmet12<-read.csv("Sotonmet12Nov2019.csv",
                   stringsAsFactors = FALSE)
sotmet13<-read.csv("Sotonmet13Nov2019.csv",
                   stringsAsFactors = FALSE)
sotmet14<-read.csv("Sot14Nov2019.csv",
                   stringsAsFactors = FALSE)
sotmet15<-read.csv("Sot15Nov2019.csv",
                   stringsAsFactors = FALSE)

#merge both sotonmet files
sotmet<- rbind(sotmet11,sotmet12,sotmet13,sotmet14,sotmet15)

# convert Date_Time into the same format for each table
UoStemp2019$Date.Time <- dmy_hm(UoStemp2019$Date.Time)
#I'm now replacing the Date.Time column head in UoStemp2019 to Date_Time
names(UoStemp2019)[2]  <- "Date_Time"
UoStemp2019 = subset(UoStemp2019, select = -c(Date,Time) )#xts doesn't work when there are Data Time columns that aren't in the right class

caltemp$Date_Time <- mdy_hms(caltemp$Date_Time)

#in this line I am also pasting the Date and time columns into a single Date_Time coloumn
sotmet$Date_Time <- dmy_hm(paste(sotmet$Date, sotmet$Time))
sotmet = subset(sotmet, select = -c(Date,Time) )

UKHO$Date_Time <- dmy_hm(UKHO$Date_Time)

#Import the data into a time series object

caltemp.ts = as.xts(caltemp, frequency = 60*24) # measurements were taken every minute and I'm setting my season as 1 day
caltemp.ts.interp = na.approx(object = caltemp.ts,
                              xout = seq(min(caltemp$Date_Time), max(caltemp$Date_Time), "min")) # using linear interpolation to fill in the missing values
plot(caltemp.ts.interp$Temperature..Deg.C., col = 2)

lines(caltemp.ts$Temperature..Deg.C., col = 1) #red section shows were the miss data has been filling in

sotmet.ts = as.xts(sotmet, frequency = 12*24) # these observations were every 5 minutes
sotmet.ts.interp = na.approx(object = sotmet.ts,
                              xout = seq(min(sotmet$Date_Time), max(sotmet$Date_Time), "5 min")) #adds approximation to missing values


UKHO.ts = as.xts(UKHO$Tide_Height.m., order.by = UKHO$Date_Time, frequency = 60*24) # these observations were every minute
plot(sotmet.ts.interp$DEPTH, col = 2)

lines(UKHO.ts, col = 1) #comparing tide hieghts data

UoStemp2019.ts = as.xts(UoStemp2019, frequency = 4*24)# these observations were every 15 minutes
plot(UoStemp2019.ts)

lines(caltemp.ts.interp$Temperature..Deg.C., col = 2) #comparing temperature data

# Converting all to 15 min intervals
timeLim<- as.POSIXct(strptime(c("2019-11-11 12:00:00", "2019-11-15 09:00:00"), format = "%Y-%m-%d %H:%M:%S"))
UKHO.15.ts <- rollmean(UKHO.ts, 15, align = "right")
plot(UKHO.ts, col=1)

lines(UKHO.15.ts, col=2) 

UKHO.15.ts<-fortify(UKHO.15.ts)
tidePlot <- ggplot(UKHO.15.ts, aes(x=Index, y=UKHO.15.ts)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab("Tide (m)") +
  theme_q2r()
tidePlot
## Warning: Removed 1605 rows containing missing values (`geom_line()`).

caltemp.i15.ts <- rollmean(caltemp.ts.interp, 15, align = "right")
plot(caltemp.i15.ts$Temperature..Deg.C., col=1)

lines(caltemp.ts.interp$Temperature..Deg.C., col=2)

caltemp.i15.ts <- fortify(caltemp.i15.ts)
tempPlot <- ggplot(caltemp.i15.ts, aes(x=Index, y=Temperature..Deg.C.)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab("Water (°C)") +
  theme_q2r()
tempPlot
## Warning: Removed 1605 rows containing missing values (`geom_line()`).

plot(UoStemp2019.ts, col=2)
lines(caltemp.i15.ts$Temperature..Deg.C., col=1)

sotATMP.15.ts <- rollmean(sotmet.ts.interp$ATMP, 3, align = "right")
plot(sotmet.ts.interp$ATMP, col=1)

lines(sotATMP.15.ts, col=2)

sotATMP.15.ts <- fortify(sotATMP.15.ts)
airTempPlot <- ggplot(sotATMP.15.ts, aes(x=Index, y=ATMP)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab("Air (°C)") +
  theme_q2r()
airTempPlot
## Warning: Removed 321 rows containing missing values (`geom_line()`).

sotWSPD.15.ts <- rollmean(sotmet.ts.interp$WSPD, 3, align = "right")
plot(sotmet.ts.interp$WSPD, col=1)

lines(sotWSPD.15.ts, col=2)

sotWSPD.15.ts <- fortify(sotWSPD.15.ts)
windPlot <- ggplot(sotWSPD.15.ts, aes(x=Index, y=WSPD)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab("Wind (knots)") +
  theme_q2r()
windPlot
## Warning: Removed 321 rows containing missing values (`geom_line()`).

sotWD.15.ts <- rollmean(sotmet.ts.interp$WD, 3, align = "right")
plot(sotmet.ts.interp$WD, col=1)

lines(sotWD.15.ts, col=2)

sotGST.15.ts <- rollmean(sotmet.ts.interp$GST, 3, align = "right")
plot(sotmet.ts.interp$GST, col=1)

lines(sotGST.15.ts, col=2)

gustPlot <- ggplot(sotGST.15.ts, aes(x=Index, y=GST)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab("Gust (knots)") +
  theme_q2r()
gustPlot
## Warning: Removed 321 rows containing missing values (`geom_line()`).

sotBARO.15.ts <- rollmean(sotmet.ts.interp$BARO, 3, align = "right")
plot(sotmet.ts.interp$BARO, col=1)

lines(sotBARO.15.ts, col=2)

sotBARO.15.ts <- fortify(sotBARO.15.ts)
presurePlot <- ggplot(sotBARO.15.ts, aes(x=Index, y=BARO)) +
  geom_line() + 
  scale_x_datetime(limits = timeLim, date_breaks = "1 day", date_labels = "%e %b") +
  xlab("Time") +
  ylab(expression(P[b]*" (mbars)")) +
  theme_q2r()
presurePlot
## Warning: Removed 321 rows containing missing values (`geom_line()`).

# Plot together

envPlot <- ggarrange(tidePlot, tempPlot, airTempPlot, presurePlot, windPlot,
                     ncol = 1,
                     align = "hv")
## Warning: Removed 1605 rows containing missing values (`geom_line()`).
## Removed 1605 rows containing missing values (`geom_line()`).
## Warning: Removed 321 rows containing missing values (`geom_line()`).
## Removed 321 rows containing missing values (`geom_line()`).
## Removed 321 rows containing missing values (`geom_line()`).
envPlot

ggsave("EnvTimeSeries.pdf", height = 5, width = 3.5, device = "pdf")

Combined Plots

Beta Diversity

nMDSPlot<- ggarrange(nMDS16S, nMDS18S,
                      nrow = 1,
                      align = "hv",
                      common.legend = TRUE,
                      legend = "bottom")
nMDSPlot

ggsave("nMDS.pdf", height=4, width=7, device="pdf")

Alpha Diversity Time Series

AlphaTime<- ggarrange(Alpha16STime, Alpha18STime,
                      ncol = 2,
                      align = "hv")
AlphaTime

ggsave("AlphaDiversity.pdf", height=4.7, width=7, device="pdf")

Stacked Taxa Time Series

speciesStacks<- ggarrange(species16Sstack, species18SStack, 
                          ncol = 1,
                          align = "hv", 
                          legend = "bottom")
speciesStacks

ggsave("Specie10TimeSeries.pdf", height=10, width=6, device="pdf")

Tables

library("sjPlot")
## Install package "strengejacke" from GitHub (`devtools::install_github("strengejacke/strengejacke")`) to load all sj-packages at once!
setwd("~/OneDrive - University of Southampton/PhD/Soton/Bioinformatics/16S")
alphaMeta<- read_tsv("alpha16S-metadata.tsv")
## Rows: 39 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## dbl (24): Sample, Time-Since-Start, UKHO_Tide_Height, Calista_Temp_oC, Soton...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
alphaEnvCorr<- tab_corr(alphaMeta, corr.method = "spearman", show.p = TRUE, p.numeric = TRUE)
alphaEnvCorr
  Sample Time-Since-Start UKHO_Tide_Height Calista_Temp_oC SotonMet_Air_Temp_oC SotonMet_Wind_Speed_knots SotonMet_Wind_Direction_degrees SotonMet_Max_Gust_knots SotonMet_Barometric_Pressure_millibars UoS_Temp_oC faith_pd pielou_evenness observed_features shannon_entropy Bacteria; Proteobacteria;
Alphaproteobacteria; SAR11_clade;
Clade_I; Clade_Ia; NA 		Bacteria; Proteobacteria;
Alphaproteobacteria; Rhodobacterales;
Rhodobacteraceae; Planktomarina; NA 		Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Cryomorphaceae;
uncultured; NA 		Bacteria; Proteobacteria;
Gammaproteobacteria; Pseudomonadales;
Thioglobaceae; SUP05_cluster; NA 		Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Flavobacteriaceae;
NS5_marine_group; NA 		Bacteria; Proteobacteria;
Gammaproteobacteria; Pseudomonadales;
Saccharospirillaceae; Oleispira; NA 		Bacteria; Proteobacteria;
Gammaproteobacteria; Burkholderiales;
Methylophilaceae; OM43_clade; NA 		Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Flavobacteriaceae;
NS3a_marine_group; uncultured_bacterium 		Bacteria; Proteobacteria;
Gammaproteobacteria; Enterobacterales;
Alteromonadaceae; Glaciecola; NA 		Bacteria; Proteobacteria;
Alphaproteobacteria; Puniceispirillales;
SAR116_clade;
Candidatus_Puniceispirillum;
uncultured_bacterium
Sample   -1.000
(<.001) 		0.000
(1.00) 		0.805
(<.001) 		0.257
(.115) 		0.198
(.227) 		0.436
(.006) 		0.207
(.207) 		0.076
(.645) 		0.883
(<.001) 		0.025
(.879) 		0.274
(.092) 		0.075
(.649) 		0.233
(.152) 		0.512
(.001) 		0.004
(.979) 		0.657
(<.001) 		0.411
(.009) 		0.306
(.058) 		-0.695
(<.001) 		0.362
(.024) 		0.655
(<.001) 		0.496
(.001) 		0.532
(<.001)
Time-Since-Start -1.000
(<.001) 		  0.000
(1.00) 		-0.805
(<.001) 		-0.257
(.115) 		-0.198
(.227) 		-0.436
(.006) 		-0.207
(.207) 		-0.076
(.645) 		-0.883
(<.001) 		-0.025
(.879) 		-0.274
(.092) 		-0.075
(.649) 		-0.233
(.152) 		-0.512
(.001) 		-0.004
(.979) 		-0.657
(<.001) 		-0.411
(.009) 		-0.306
(.058) 		0.695
(<.001) 		-0.362
(.024) 		-0.655
(<.001) 		-0.496
(.001) 		-0.532
(<.001)
UKHO_Tide_Height 0.000
(1.00) 		0.000
(1.00) 		  0.204
(.212) 		0.024
(.887) 		0.106
(.520) 		0.228
(.162) 		0.052
(.754) 		0.021
(.901) 		-0.198
(.227) 		-0.150
(.361) 		-0.221
(.177) 		-0.072
(.662) 		-0.219
(.180) 		-0.018
(.911) 		-0.109
(.509) 		-0.290
(.073) 		0.226
(.167) 		0.062
(.709) 		0.031
(.853) 		0.184
(.262) 		-0.300
(.064) 		-0.161
(.326) 		0.112
(.496)
Calista_Temp_oC 0.805
(<.001) 		-0.805
(<.001) 		0.204
(.212) 		  0.426
(.007) 		0.116
(.482) 		0.487
(.002) 		0.118
(.473) 		0.057
(.730) 		0.727
(<.001) 		-0.155
(.343) 		0.280
(.085) 		-0.132
(.423) 		0.008
(.963) 		0.371
(.021) 		-0.013
(.939) 		0.456
(.004) 		0.244
(.135) 		0.187
(.254) 		-0.486
(.002) 		0.177
(.281) 		0.472
(.003) 		0.272
(.094) 		0.338
(.035)
SotonMet_Air_Temp_oC 0.257
(.115) 		-0.257
(.115) 		0.024
(.887) 		0.426
(.007) 		  0.368
(.021) 		0.227
(.165) 		0.386
(.015) 		-0.216
(.187) 		0.442
(.005) 		0.015
(.926) 		0.262
(.108) 		-0.172
(.295) 		0.023
(.888) 		0.223
(.172) 		-0.041
(.805) 		0.091
(.582) 		0.059
(.720) 		-0.129
(.435) 		-0.125
(.447) 		-0.027
(.872) 		0.010
(.953) 		-0.081
(.624) 		-0.146
(.376)
SotonMet_Wind_Speed_knots 0.198
(.227) 		-0.198
(.227) 		0.106
(.520) 		0.116
(.482) 		0.368
(.021) 		  -0.022
(.895) 		0.961
(<.001) 		0.217
(.185) 		0.260
(.110) 		0.094
(.570) 		-0.352
(.028) 		0.017
(.917) 		-0.132
(.423) 		-0.069
(.674) 		0.367
(.021) 		0.253
(.120) 		0.386
(.015) 		-0.130
(.429) 		-0.348
(.030) 		-0.054
(.743) 		0.120
(.468) 		0.279
(.085) 		0.131
(.428)
SotonMet_Wind_Direction_degrees 0.436
(.006) 		-0.436
(.006) 		0.228
(.162) 		0.487
(.002) 		0.227
(.165) 		-0.022
(.895) 		  0.030
(.857) 		0.060
(.717) 		0.296
(.067) 		0.080
(.627) 		-0.002
(.989) 		-0.025
(.881) 		0.019
(.911) 		0.432
(.006) 		-0.049
(.769) 		0.217
(.186) 		0.260
(.110) 		0.305
(.059) 		-0.471
(.002) 		0.292
(.071) 		0.157
(.341) 		0.033
(.841) 		0.347
(.031)
SotonMet_Max_Gust_knots 0.207
(.207) 		-0.207
(.207) 		0.052
(.754) 		0.118
(.473) 		0.386
(.015) 		0.961
(<.001) 		0.030
(.857) 		  0.218
(.182) 		0.292
(.071) 		0.100
(.546) 		-0.318
(.049) 		-0.045
(.785) 		-0.152
(.355) 		-0.036
(.829) 		0.389
(.014) 		0.294
(.069) 		0.372
(.020) 		-0.064
(.701) 		-0.367
(.022) 		-0.062
(.708) 		0.191
(.244) 		0.359
(.025) 		0.130
(.429)
SotonMet_Barometric_Pressure_millibars 0.076
(.645) 		-0.076
(.645) 		0.021
(.901) 		0.057
(.730) 		-0.216
(.187) 		0.217
(.185) 		0.060
(.717) 		0.218
(.182) 		  -0.060
(.715) 		-0.459
(.003) 		-0.182
(.267) 		-0.318
(.048) 		-0.357
(.026) 		0.054
(.742) 		0.456
(.003) 		0.282
(.082) 		0.437
(.005) 		0.201
(.220) 		-0.495
(.001) 		0.038
(.818) 		0.150
(.363) 		0.333
(.038) 		0.274
(.091)
UoS_Temp_oC 0.883
(<.001) 		-0.883
(<.001) 		-0.198
(.227) 		0.727
(<.001) 		0.442
(.005) 		0.260
(.110) 		0.296
(.067) 		0.292
(.071) 		-0.060
(.715) 		  0.089
(.589) 		0.301
(.063) 		0.018
(.915) 		0.195
(.235) 		0.389
(.014) 		0.064
(.697) 		0.598
(<.001) 		0.231
(.157) 		0.124
(.450) 		-0.582
(<.001) 		0.148
(.367) 		0.586
(<.001) 		0.441
(.005) 		0.323
(.045)
faith_pd 0.025
(.879) 		-0.025
(.879) 		-0.150
(.361) 		-0.155
(.343) 		0.015
(.926) 		0.094
(.570) 		0.080
(.627) 		0.100
(.546) 		-0.459
(.003) 		0.089
(.589) 		  -0.125
(.447) 		0.572
(<.001) 		0.474
(.003) 		-0.032
(.845) 		-0.121
(.461) 		-0.101
(.541) 		-0.235
(.150) 		-0.072
(.662) 		0.217
(.185) 		-0.040
(.806) 		0.017
(.916) 		-0.125
(.448) 		0.043
(.795)
pielou_evenness 0.274
(.092) 		-0.274
(.092) 		-0.221
(.177) 		0.280
(.085) 		0.262
(.108) 		-0.352
(.028) 		-0.002
(.989) 		-0.318
(.049) 		-0.182
(.267) 		0.301
(.063) 		-0.125
(.447) 		  -0.166
(.312) 		0.377
(.019) 		0.375
(.019) 		-0.229
(.161) 		0.178
(.277) 		-0.040
(.811) 		0.174
(.288) 		-0.023
(.890) 		0.238
(.144) 		0.277
(.088) 		0.051
(.756) 		0.132
(.424)
observed_features 0.075
(.649) 		-0.075
(.649) 		-0.072
(.662) 		-0.132
(.423) 		-0.172
(.295) 		0.017
(.917) 		-0.025
(.881) 		-0.045
(.785) 		-0.318
(.048) 		0.018
(.915) 		0.572
(<.001) 		-0.166
(.312) 		  0.804
(<.001) 		-0.019
(.910) 		-0.126
(.444) 		-0.023
(.890) 		0.082
(.618) 		-0.082
(.621) 		0.059
(.723) 		-0.045
(.787) 		-0.022
(.895) 		-0.018
(.913) 		0.196
(.232)
shannon_entropy 0.233
(.152) 		-0.233
(.152) 		-0.219
(.180) 		0.008
(.963) 		0.023
(.888) 		-0.132
(.423) 		0.019
(.911) 		-0.152
(.355) 		-0.357
(.026) 		0.195
(.235) 		0.474
(.003) 		0.377
(.019) 		0.804
(<.001) 		  0.198
(.226) 		-0.186
(.256) 		0.135
(.410) 		0.132
(.423) 		0.021
(.901) 		-0.022
(.893) 		0.119
(.471) 		0.152
(.355) 		0.059
(.723) 		0.273
(.093)
Bacteria; Proteobacteria;
Alphaproteobacteria; SAR11_clade;
Clade_I; Clade_Ia; NA 		0.512
(.001) 		-0.512
(.001) 		-0.018
(.911) 		0.371
(.021) 		0.223
(.172) 		-0.069
(.674) 		0.432
(.006) 		-0.036
(.829) 		0.054
(.742) 		0.389
(.014) 		-0.032
(.845) 		0.375
(.019) 		-0.019
(.910) 		0.198
(.226) 		  -0.493
(.001) 		0.095
(.564) 		0.426
(.007) 		0.693
(<.001) 		-0.390
(.014) 		0.820
(<.001) 		0.146
(.374) 		0.017
(.919) 		0.481
(.002)
Bacteria; Proteobacteria;
Alphaproteobacteria; Rhodobacterales;
Rhodobacteraceae; Planktomarina; NA 		0.004
(.979) 		-0.004
(.979) 		-0.109
(.509) 		-0.013
(.939) 		-0.041
(.805) 		0.367
(.021) 		-0.049
(.769) 		0.389
(.014) 		0.456
(.003) 		0.064
(.697) 		-0.121
(.461) 		-0.229
(.161) 		-0.126
(.444) 		-0.186
(.256) 		-0.493
(.001) 		  0.351
(.028) 		0.144
(.381) 		-0.290
(.073) 		-0.317
(.049) 		-0.517
(.001) 		0.214
(.191) 		0.332
(.039) 		0.044
(.792)
Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Cryomorphaceae;
uncultured; NA 		0.657
(<.001) 		-0.657
(<.001) 		-0.290
(.073) 		0.456
(.004) 		0.091
(.582) 		0.253
(.120) 		0.217
(.186) 		0.294
(.069) 		0.282
(.082) 		0.598
(<.001) 		-0.101
(.541) 		0.178
(.277) 		-0.023
(.890) 		0.135
(.410) 		0.095
(.564) 		0.351
(.028) 		  0.396
(.013) 		0.056
(.735) 		-0.755
(<.001) 		-0.007
(.967) 		0.903
(<.001) 		0.817
(<.001) 		0.392
(.014)
Bacteria; Proteobacteria;
Gammaproteobacteria; Pseudomonadales;
Thioglobaceae; SUP05_cluster; NA 		0.411
(.009) 		-0.411
(.009) 		0.226
(.167) 		0.244
(.135) 		0.059
(.720) 		0.386
(.015) 		0.260
(.110) 		0.372
(.020) 		0.437
(.005) 		0.231
(.157) 		-0.235
(.150) 		-0.040
(.811) 		0.082
(.618) 		0.132
(.423) 		0.426
(.007) 		0.144
(.381) 		0.396
(.013) 		  0.284
(.080) 		-0.689
(<.001) 		0.497
(.001) 		0.228
(.162) 		0.468
(.003) 		0.685
(<.001)
Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Flavobacteriaceae;
NS5_marine_group; NA 		0.306
(.058) 		-0.306
(.058) 		0.062
(.709) 		0.187
(.254) 		-0.129
(.435) 		-0.130
(.429) 		0.305
(.059) 		-0.064
(.701) 		0.201
(.220) 		0.124
(.450) 		-0.072
(.662) 		0.174
(.288) 		-0.082
(.621) 		0.021
(.901) 		0.693
(<.001) 		-0.290
(.073) 		0.056
(.735) 		0.284
(.080) 		  -0.348
(.030) 		0.740
(<.001) 		0.243
(.137) 		0.170
(.301) 		0.565
(<.001)
Bacteria; Proteobacteria;
Gammaproteobacteria; Pseudomonadales;
Saccharospirillaceae; Oleispira; NA 		-0.695
(<.001) 		0.695
(<.001) 		0.031
(.853) 		-0.486
(.002) 		-0.125
(.447) 		-0.348
(.030) 		-0.471
(.002) 		-0.367
(.022) 		-0.495
(.001) 		-0.582
(<.001) 		0.217
(.185) 		-0.023
(.890) 		0.059
(.723) 		-0.022
(.893) 		-0.390
(.014) 		-0.317
(.049) 		-0.755
(<.001) 		-0.689
(<.001) 		-0.348
(.030) 		  -0.318
(.049) 		-0.604
(<.001) 		-0.667
(<.001) 		-0.536
(<.001)
Bacteria; Proteobacteria;
Gammaproteobacteria; Burkholderiales;
Methylophilaceae; OM43_clade; NA 		0.362
(.024) 		-0.362
(.024) 		0.184
(.262) 		0.177
(.281) 		-0.027
(.872) 		-0.054
(.743) 		0.292
(.071) 		-0.062
(.708) 		0.038
(.818) 		0.148
(.367) 		-0.040
(.806) 		0.238
(.144) 		-0.045
(.787) 		0.119
(.471) 		0.820
(<.001) 		-0.517
(.001) 		-0.007
(.967) 		0.497
(.001) 		0.740
(<.001) 		-0.318
(.049) 		  0.045
(.783) 		-0.001
(.993) 		0.494
(.001)
Bacteria; Bacteroidota; Bacteroidia;
Flavobacteriales; Flavobacteriaceae;
NS3a_marine_group; uncultured_bacterium 		0.655
(<.001) 		-0.655
(<.001) 		-0.300
(.064) 		0.472
(.003) 		0.010
(.953) 		0.120
(.468) 		0.157
(.341) 		0.191
(.244) 		0.150
(.363) 		0.586
(<.001) 		0.017
(.916) 		0.277
(.088) 		-0.022
(.895) 		0.152
(.355) 		0.146
(.374) 		0.214
(.191) 		0.903
(<.001) 		0.228
(.162) 		0.243
(.137) 		-0.604
(<.001) 		0.045
(.783) 		  0.793
(<.001) 		0.457
(.003)
Bacteria; Proteobacteria;
Gammaproteobacteria; Enterobacterales;
Alteromonadaceae; Glaciecola; NA 		0.496
(.001) 		-0.496
(.001) 		-0.161
(.326) 		0.272
(.094) 		-0.081
(.624) 		0.279
(.085) 		0.033
(.841) 		0.359
(.025) 		0.333
(.038) 		0.441
(.005) 		-0.125
(.448) 		0.051
(.756) 		-0.018
(.913) 		0.059
(.723) 		0.017
(.919) 		0.332
(.039) 		0.817
(<.001) 		0.468
(.003) 		0.170
(.301) 		-0.667
(<.001) 		-0.001
(.993) 		0.793
(<.001) 		  0.438
(.005)
Bacteria; Proteobacteria;
Alphaproteobacteria; Puniceispirillales;
SAR116_clade;
Candidatus_Puniceispirillum;
uncultured_bacterium 		0.532
(<.001) 		-0.532
(<.001) 		0.112
(.496) 		0.338
(.035) 		-0.146
(.376) 		0.131
(.428) 		0.347
(.031) 		0.130
(.429) 		0.274
(.091) 		0.323
(.045) 		0.043
(.795) 		0.132
(.424) 		0.196
(.232) 		0.273
(.093) 		0.481
(.002) 		0.044
(.792) 		0.392
(.014) 		0.685
(<.001) 		0.565
(<.001) 		-0.536
(<.001) 		0.494
(.001) 		0.457
(.003) 		0.438
(.005) 		 
Computed correlation used spearman-method with listwise-deletion. 
setwd("~/OneDrive - University of Southampton/PhD/Soton/Bioinformatics/18S_dada2redo")
alphaMeta<- read_tsv("alpha18S-metadata.tsv")
## Rows: 39 Columns: 24
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: "\t"
## dbl (24): id, Time-Since-Start, UKHO_Tide_Height, Calista_Temp_oC, SotonMet_...
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
alphaEnvCorr<- tab_corr(alphaMeta, corr.method = "spearman", show.p = TRUE, p.numeric = TRUE)
alphaEnvCorr
  id Time-Since-Start UKHO_Tide_Height Calista_Temp_oC SotonMet_Air_Temp_oC SotonMet_Wind_Speed_knots SotonMet_Wind_Direction_degrees SotonMet_Max_Gust_knots SotonMet_Barometric_Pressure_millibars UoS_Temp_oC faith_pd pielou_evenness observed_features shannon_entropy Eukaryota; NA; NA; NA; NA; NA; NA Eukaryota; Arthropoda; Maxillopoda; NA;
NA; NA; NA 		Eukaryota; Annelida; Polychaeta;
Capitellida; Capitellida; Capitellida;
NA 		Eukaryota; Dinoflagellata; Dinophyceae;
NA; NA; NA; NA 		Eukaryota; Cryptophyceae; Cryptophyceae;
Cryptomonadales; Cryptomonadales;
uncultured; Cryptophyta_sp. 		Eukaryota; Picozoa; Picomonadida;
Picomonadida; Picomonadida;
Picomonadida; uncultured_phototrophic 		Eukaryota; Diatomea; Mediophyceae;
Mediophyceae; Mediophyceae;
Thalassiosira; NA 		Eukaryota; Protalveolata; Syndiniales;
Syndiniales; Syndiniales_Group_I;
Syndiniales_Group_I;
Karlodinium_veneficum 		Eukaryota; Arthropoda; Maxillopoda;
Calanoida; Calanoida; Calanoida;
uncultured_eukaryote 		Eukaryota; Incertae_Sedis;
Incertae_Sedis; Incertae_Sedis;
Incertae_Sedis; Telonema;
Telonema_antarcticum
id   -1.000
(<.001) 		0.000
(1.00) 		0.805
(<.001) 		0.257
(.115) 		0.198
(.227) 		0.436
(.006) 		0.207
(.207) 		0.076
(.645) 		0.883
(<.001) 		-0.048
(.772) 		-0.305
(.060) 		-0.349
(.030) 		-0.278
(.087) 		0.050
(.764) 		-0.168
(.305) 		0.314
(.052) 		0.195
(.233) 		-0.309
(.056) 		0.028
(.864) 		-0.209
(.201) 		-0.225
(.168) 		0.339
(.035) 		-0.031
(.852)
Time-Since-Start -1.000
(<.001) 		  0.000
(1.00) 		-0.805
(<.001) 		-0.257
(.115) 		-0.198
(.227) 		-0.436
(.006) 		-0.207
(.207) 		-0.076
(.645) 		-0.883
(<.001) 		0.048
(.772) 		0.305
(.060) 		0.349
(.030) 		0.278
(.087) 		-0.050
(.764) 		0.168
(.305) 		-0.314
(.052) 		-0.195
(.233) 		0.309
(.056) 		-0.028
(.864) 		0.209
(.201) 		0.225
(.168) 		-0.339
(.035) 		0.031
(.852)
UKHO_Tide_Height 0.000
(1.00) 		0.000
(1.00) 		  0.204
(.212) 		0.024
(.887) 		0.106
(.520) 		0.228
(.162) 		0.052
(.754) 		0.021
(.901) 		-0.198
(.227) 		-0.001
(.995) 		0.046
(.783) 		-0.079
(.635) 		-0.020
(.903) 		0.140
(.394) 		0.185
(.259) 		-0.078
(.638) 		-0.088
(.593) 		0.025
(.880) 		0.069
(.677) 		0.215
(.190) 		0.113
(.491) 		0.372
(.020) 		0.238
(.144)
Calista_Temp_oC 0.805
(<.001) 		-0.805
(<.001) 		0.204
(.212) 		  0.426
(.007) 		0.116
(.482) 		0.487
(.002) 		0.118
(.473) 		0.057
(.730) 		0.727
(<.001) 		-0.325
(.044) 		-0.246
(.131) 		-0.544
(<.001) 		-0.327
(.043) 		0.159
(.331) 		-0.155
(.343) 		0.353
(.028) 		0.115
(.483) 		-0.455
(.004) 		0.111
(.500) 		-0.121
(.462) 		-0.008
(.962) 		0.520
(.001) 		0.053
(.749)
SotonMet_Air_Temp_oC 0.257
(.115) 		-0.257
(.115) 		0.024
(.887) 		0.426
(.007) 		  0.368
(.021) 		0.227
(.165) 		0.386
(.015) 		-0.216
(.187) 		0.442
(.005) 		-0.151
(.360) 		-0.151
(.357) 		-0.361
(.024) 		-0.249
(.127) 		0.118
(.474) 		0.046
(.782) 		0.133
(.420) 		0.080
(.628) 		-0.288
(.075) 		0.030
(.857) 		-0.093
(.572) 		0.170
(.301) 		0.262
(.107) 		-0.039
(.814)
SotonMet_Wind_Speed_knots 0.198
(.227) 		-0.198
(.227) 		0.106
(.520) 		0.116
(.482) 		0.368
(.021) 		  -0.022
(.895) 		0.961
(<.001) 		0.217
(.185) 		0.260
(.110) 		0.097
(.557) 		-0.202
(.218) 		-0.166
(.312) 		-0.220
(.179) 		0.021
(.897) 		-0.001
(.997) 		0.028
(.866) 		0.078
(.638) 		-0.043
(.796) 		0.036
(.828) 		-0.107
(.518) 		-0.220
(.178) 		-0.082
(.622) 		-0.144
(.383)
SotonMet_Wind_Direction_degrees 0.436
(.006) 		-0.436
(.006) 		0.228
(.162) 		0.487
(.002) 		0.227
(.165) 		-0.022
(.895) 		  0.030
(.857) 		0.060
(.717) 		0.296
(.067) 		-0.122
(.459) 		-0.012
(.941) 		-0.298
(.065) 		-0.113
(.492) 		0.041
(.804) 		0.042
(.799) 		0.027
(.869) 		0.403
(.011) 		-0.014
(.934) 		-0.031
(.850) 		-0.190
(.247) 		-0.079
(.632) 		0.475
(.002) 		-0.039
(.815)
SotonMet_Max_Gust_knots 0.207
(.207) 		-0.207
(.207) 		0.052
(.754) 		0.118
(.473) 		0.386
(.015) 		0.961
(<.001) 		0.030
(.857) 		  0.218
(.182) 		0.292
(.071) 		0.177
(.280) 		-0.224
(.171) 		-0.198
(.227) 		-0.237
(.146) 		0.042
(.800) 		-0.014
(.930) 		0.079
(.632) 		0.106
(.520) 		-0.054
(.746) 		0.045
(.784) 		-0.110
(.504) 		-0.187
(.255) 		-0.045
(.787) 		-0.120
(.467)
SotonMet_Barometric_Pressure_millibars 0.076
(.645) 		-0.076
(.645) 		0.021
(.901) 		0.057
(.730) 		-0.216
(.187) 		0.217
(.185) 		0.060
(.717) 		0.218
(.182) 		  -0.060
(.715) 		-0.232
(.155) 		0.157
(.341) 		-0.220
(.178) 		0.108
(.513) 		-0.247
(.130) 		0.033
(.842) 		-0.121
(.464) 		0.294
(.069) 		0.041
(.803) 		-0.095
(.564) 		-0.132
(.424) 		-0.313
(.053) 		0.079
(.632) 		-0.105
(.525)
UoS_Temp_oC 0.883
(<.001) 		-0.883
(<.001) 		-0.198
(.227) 		0.727
(<.001) 		0.442
(.005) 		0.260
(.110) 		0.296
(.067) 		0.292
(.071) 		-0.060
(.715) 		  -0.078
(.638) 		-0.319
(.047) 		-0.402
(.011) 		-0.312
(.053) 		-0.013
(.937) 		-0.094
(.571) 		0.310
(.055) 		0.049
(.767) 		-0.361
(.024) 		-0.093
(.575) 		-0.209
(.202) 		-0.256
(.116) 		0.283
(.080) 		-0.160
(.330)
faith_pd -0.048
(.772) 		0.048
(.772) 		-0.001
(.995) 		-0.325
(.044) 		-0.151
(.360) 		0.097
(.557) 		-0.122
(.459) 		0.177
(.280) 		-0.232
(.155) 		-0.078
(.638) 		  0.202
(.216) 		0.678
(<.001) 		0.386
(.016) 		-0.082
(.618) 		-0.086
(.600) 		-0.110
(.502) 		0.081
(.624) 		0.385
(.016) 		-0.050
(.761) 		0.027
(.871) 		-0.054
(.744) 		-0.288
(.076) 		0.101
(.541)
pielou_evenness -0.305
(.060) 		0.305
(.060) 		0.046
(.783) 		-0.246
(.131) 		-0.151
(.357) 		-0.202
(.218) 		-0.012
(.941) 		-0.224
(.171) 		0.157
(.341) 		-0.319
(.047) 		0.202
(.216) 		  0.381
(.017) 		0.934
(<.001) 		-0.073
(.659) 		-0.011
(.949) 		-0.198
(.227) 		0.171
(.296) 		0.391
(.014) 		-0.090
(.585) 		0.050
(.764) 		0.038
(.816) 		-0.080
(.629) 		0.100
(.544)
observed_features -0.349
(.030) 		0.349
(.030) 		-0.079
(.635) 		-0.544
(<.001) 		-0.361
(.024) 		-0.166
(.312) 		-0.298
(.065) 		-0.198
(.227) 		-0.220
(.178) 		-0.402
(.011) 		0.678
(<.001) 		0.381
(.017) 		  0.634
(<.001) 		-0.087
(.598) 		-0.091
(.583) 		-0.341
(.034) 		-0.066
(.688) 		0.452
(.004) 		-0.118
(.476) 		0.000
(1.00) 		-0.043
(.797) 		-0.473
(.002) 		-0.009
(.957)
shannon_entropy -0.278
(.087) 		0.278
(.087) 		-0.020
(.903) 		-0.327
(.043) 		-0.249
(.127) 		-0.220
(.179) 		-0.113
(.492) 		-0.237
(.146) 		0.108
(.513) 		-0.312
(.053) 		0.386
(.016) 		0.934
(<.001) 		0.634
(<.001) 		  -0.092
(.577) 		-0.067
(.684) 		-0.209
(.201) 		0.125
(.448) 		0.438
(.005) 		-0.105
(.525) 		0.038
(.819) 		0.002
(.990) 		-0.221
(.176) 		0.085
(.609)
Eukaryota; NA; NA; NA; NA; NA; NA 0.050
(.764) 		-0.050
(.764) 		0.140
(.394) 		0.159
(.331) 		0.118
(.474) 		0.021
(.897) 		0.041
(.804) 		0.042
(.800) 		-0.247
(.130) 		-0.013
(.937) 		-0.082
(.618) 		-0.073
(.659) 		-0.087
(.598) 		-0.092
(.577) 		  -0.456
(.004) 		0.016
(.925) 		0.056
(.733) 		-0.427
(.007) 		0.846
(<.001) 		0.455
(.004) 		0.789
(<.001) 		0.233
(.153) 		0.668
(<.001)
Eukaryota; Arthropoda; Maxillopoda; NA;
NA; NA; NA 		-0.168
(.305) 		0.168
(.305) 		0.185
(.259) 		-0.155
(.343) 		0.046
(.782) 		-0.001
(.997) 		0.042
(.799) 		-0.014
(.930) 		0.033
(.842) 		-0.094
(.571) 		-0.086
(.600) 		-0.011
(.949) 		-0.091
(.583) 		-0.067
(.684) 		-0.456
(.004) 		  -0.178
(.277) 		-0.095
(.563) 		0.011
(.945) 		-0.522
(.001) 		-0.199
(.224) 		-0.323
(.045) 		-0.150
(.361) 		-0.226
(.167)
Eukaryota; Annelida; Polychaeta;
Capitellida; Capitellida; Capitellida;
NA 		0.314
(.052) 		-0.314
(.052) 		-0.078
(.638) 		0.353
(.028) 		0.133
(.420) 		0.028
(.866) 		0.027
(.869) 		0.079
(.632) 		-0.121
(.464) 		0.310
(.055) 		-0.110
(.502) 		-0.198
(.227) 		-0.341
(.034) 		-0.209
(.201) 		0.016
(.925) 		-0.178
(.277) 		  0.016
(.925) 		-0.040
(.808) 		0.123
(.456) 		0.012
(.940) 		0.011
(.947) 		0.149
(.366) 		-0.088
(.594)
Eukaryota; Dinoflagellata; Dinophyceae;
NA; NA; NA; NA 		0.195
(.233) 		-0.195
(.233) 		-0.088
(.593) 		0.115
(.483) 		0.080
(.628) 		0.078
(.638) 		0.403
(.011) 		0.106
(.520) 		0.294
(.069) 		0.049
(.767) 		0.081
(.624) 		0.171
(.296) 		-0.066
(.688) 		0.125
(.448) 		0.056
(.733) 		-0.095
(.563) 		0.016
(.925) 		  0.340
(.034) 		0.069
(.676) 		-0.229
(.162) 		0.017
(.918) 		0.003
(.986) 		0.140
(.394)
Eukaryota; Cryptophyceae; Cryptophyceae;
Cryptomonadales; Cryptomonadales;
uncultured; Cryptophyta_sp. 		-0.309
(.056) 		0.309
(.056) 		0.025
(.880) 		-0.455
(.004) 		-0.288
(.075) 		-0.043
(.796) 		-0.014
(.934) 		-0.054
(.746) 		0.041
(.803) 		-0.361
(.024) 		0.385
(.016) 		0.391
(.014) 		0.452
(.004) 		0.438
(.005) 		-0.427
(.007) 		0.011
(.945) 		-0.040
(.808) 		0.340
(.034) 		  -0.398
(.012) 		-0.049
(.768) 		-0.353
(.028) 		-0.386
(.015) 		-0.217
(.185)
Eukaryota; Picozoa; Picomonadida;
Picomonadida; Picomonadida;
Picomonadida; uncultured_phototrophic 		0.028
(.864) 		-0.028
(.864) 		0.069
(.677) 		0.111
(.500) 		0.030
(.857) 		0.036
(.828) 		-0.031
(.850) 		0.045
(.784) 		-0.095
(.564) 		-0.093
(.575) 		-0.050
(.761) 		-0.090
(.585) 		-0.118
(.476) 		-0.105
(.525) 		0.846
(<.001) 		-0.522
(.001) 		0.123
(.456) 		0.069
(.676) 		-0.398
(.012) 		  0.488
(.002) 		0.810
(<.001) 		0.110
(.505) 		0.778
(<.001)
Eukaryota; Diatomea; Mediophyceae;
Mediophyceae; Mediophyceae;
Thalassiosira; NA 		-0.209
(.201) 		0.209
(.201) 		0.215
(.190) 		-0.121
(.462) 		-0.093
(.572) 		-0.107
(.518) 		-0.190
(.247) 		-0.110
(.504) 		-0.132
(.424) 		-0.209
(.202) 		0.027
(.871) 		0.050
(.764) 		0.000
(1.00) 		0.038
(.819) 		0.455
(.004) 		-0.199
(.224) 		0.012
(.940) 		-0.229
(.162) 		-0.049
(.768) 		0.488
(.002) 		  0.539
(<.001) 		-0.118
(.476) 		0.438
(.005)
Eukaryota; Protalveolata; Syndiniales;
Syndiniales; Syndiniales_Group_I;
Syndiniales_Group_I;
Karlodinium_veneficum 		-0.225
(.168) 		0.225
(.168) 		0.113
(.491) 		-0.008
(.962) 		0.170
(.301) 		-0.220
(.178) 		-0.079
(.632) 		-0.187
(.255) 		-0.313
(.053) 		-0.256
(.116) 		-0.054
(.744) 		0.038
(.816) 		-0.043
(.797) 		0.002
(.990) 		0.789
(<.001) 		-0.323
(.045) 		0.011
(.947) 		0.017
(.918) 		-0.353
(.028) 		0.810
(<.001) 		0.539
(<.001) 		  0.117
(.477) 		0.756
(<.001)
Eukaryota; Arthropoda; Maxillopoda;
Calanoida; Calanoida; Calanoida;
uncultured_eukaryote 		0.339
(.035) 		-0.339
(.035) 		0.372
(.020) 		0.520
(.001) 		0.262
(.107) 		-0.082
(.622) 		0.475
(.002) 		-0.045
(.787) 		0.079
(.632) 		0.283
(.080) 		-0.288
(.076) 		-0.080
(.629) 		-0.473
(.002) 		-0.221
(.176) 		0.233
(.153) 		-0.150
(.361) 		0.149
(.366) 		0.003
(.986) 		-0.386
(.015) 		0.110
(.505) 		-0.118
(.476) 		0.117
(.477) 		  0.009
(.958)
Eukaryota; Incertae_Sedis;
Incertae_Sedis; Incertae_Sedis;
Incertae_Sedis; Telonema;
Telonema_antarcticum 		-0.031
(.852) 		0.031
(.852) 		0.238
(.144) 		0.053
(.749) 		-0.039
(.814) 		-0.144
(.383) 		-0.039
(.815) 		-0.120
(.467) 		-0.105
(.525) 		-0.160
(.330) 		0.101
(.541) 		0.100
(.544) 		-0.009
(.957) 		0.085
(.609) 		0.668
(<.001) 		-0.226
(.167) 		-0.088
(.594) 		0.140
(.394) 		-0.217
(.185) 		0.778
(<.001) 		0.438
(.005) 		0.756
(<.001) 		0.009
(.958) 		 
Computed correlation used spearman-method with listwise-deletion.